\chapter{R语言基础}
\section{R 语言初识}
\subsection{R语言的强大}
\begin{itemize}
	\item 各种回归：OLS、时间序列阈值模型(\texttt{tsDyn})、面板阈值(\texttt{pdR::ptm})、机器学习等
	\item 漂亮的图形：\texttt{ggplot2}\\ \url{http://www.cellyse.com/how_to_use_gggplot2_part1/}	
	\item 最新的算法：往往一篇文章就伴随着R代码的发布
	\item 输出到latex的各种方便：\texttt{stargazer}
\end{itemize}



\subsection{R 语言初识}

\begin{itemize}
	\item 如何获取帮助？\verb|RSiteSearch,?lm|
	\item 工作空间是什么？得到和设置工作空间的命令是什么？
	\item 列出当前工作空间的命令是什么？如何删除一个或多个变量？
	\item 如何保存当前工作空间的一个变量或者多个多个变量？保存后的文件后缀名是什么？
	\item 如何读取工作空间数据到当前会话中？
	\item 包是什么？如何安装包？如何载入包？
\end{itemize}


\section{创建数据}
\subsection{向量、矩阵和数组}

\begin{itemize}
	\item 向量：\lstinline|a <- c(1,2,3,4)|
	\item 矩阵：\lstinline|y <- matrix(1:20,nrow=5,ncol=4)|
	\item 数组：\texttt{?array}，较少使用。\verb|array(1:3,c(2,4,2))|
\end{itemize}

\textbf{约束：同一向量、矩阵和数组不可混杂不同模式的数据。}
\begin{itemize}
	\item 数值型(\texttt{numeric})：浮点型、整型（\texttt{int}）
	\item 字符型：\texttt{character}
	\item 逻辑型：\texttt{logic}
\end{itemize}	

\subsection{数据框、因子和列表}
	\begin{itemize}	
		\item 数据框：\texttt{?data.frame}。R语言的战斗机。
		\item \hypertarget{2-1}{因子}：\texttt{factor}，很少直接创造，往往是读入数据的时候形成。
	\end{itemize}
		\begin{verbatim}
		data("Affairs",package = 'AER') # `婚外情数据'
		Affairs$ynaffair[Affairs$affairs >0] <- 1
		Affairs$ynaffair[Affairs$affairs == 0] <- 0
		Affairs$ynaffair <- factor(Affairs$ynaffair,levels = c(0,1),labels = c('No','Yes'))		
		\end{verbatim}
			\begin{itemize}	
		\item 列表：\texttt{list}
	\end{itemize}

\subsection{创建数据：外部导入}
\begin{itemize}
	\item 如何输入excel、stata的数据？
	\begin{itemize}
		\item \texttt{xlsx}，之前要安装\texttt{jdk}。注意编码问题\texttt{UTF-8}。
		\item \texttt{xlsx}包导入数据慢，而且数据量大约超过2-3W行就经常内存溢出，另有一个包\texttt{openxlsx}，导入速度快，但是无法导入\texttt{.xls}格式数据。
		\item \texttt{foreign}：stata, SPSS。
	\end{itemize}
	\item 保存数据的两个常用命令：\lstinline|save,save.image|，使用帮助查看；
\end{itemize}



\subsection{一些常用的基础命令}
\includegraphics[scale=0.6]{2-1.png}

\subsection{字符串处理}
\verb|stringr|

\begin{verbatim}
fruit <- c("apple", "banana", "pear", "pinapple")
str_detect(fruit, "a") # 字符串有a就行
str_detect(fruit, "^a") # 开头有a
str_detect(fruit, "a$") # 截尾有a
str_detect(fruit, "[bl]") # 字符串有b或者l
\end{verbatim}

正则表达：\href{http://www.runoob.com/regexp/regexp-syntax.html}{正则表达式语法}（http://www.runoob.com/regexp/regexp-syntax.html）
\section{基本数据管理}
\subsection{常用操作}
\begin{itemize}
	\item 如何创建一个新的变量？如何在数据框中生成一个新的列？
	\item 如何重命名数据框的列名？
	\item 如何选取数据框中的指定行或指定列？
	\begin{itemize}
		\item 通过下标选：\lstinline|a[1,2],a[,1],a[,1],a[1:10,1],a[1,2:3]|
		\item 通过列名选：
		
		\lstinline|library(wooldrige);data(jtrain);|\\
		\lstinline|jtrain[,'year'];jtrain[,'fcode'];jtrain[,c('year','fcode')]|		
		\item 通过美元符号提取列：\lstinline|jtrain$year|
	\end{itemize}	
	\item 逻辑运算选子集；\\
	\texttt{library(wooldrige); data(jtrain)\\ jtrain[jtrain\$year==1987,]\\
		jtrain[jtrain\$year==1987 \& jtrain\$fcode==410032,]}
	\item 缺失值：
	\end{itemize}
\begin{lstlisting}[language=R]
NA
is.na(a)
sum(c(1,2,NA,4),na.rm=T)
\end{lstlisting}
	\begin{itemize}
\item 数据排序命令？

\lstinline|order(c(4,5,20,14,9)) \# for subscript|

\lstinline|sort(c(4,5,20,14,9)) \# for original data|
\item 数值、字符串、矩阵，数据框之间的转换\\
\lstinline|as.numeric; as.matrix; as.data.frame|
\end{itemize}


\subsection{高级操作}
\subsubsection{数据合并: }
	\begin{itemize}
		\item \texttt{example(merge)}
		\item \texttt{reshape2}包：\texttt{melt, dcast}，见\verb|MeltCast.R|
		
		\includegraphics[scale=0.5]{BaseofR1.png}
	\end{itemize}	

\subsubsection{时间处理}

\paragraph{日期与字符串之间的转换；}
\begin{verbatim}
strDate <- as.Date(c('2018-1-20','2018-1-22')) # 一定要包含月和日
as.character(strDate)
difftime(Sys.Date(),as.Date('2018-1-19'),units='weeks')
\end{verbatim}
\paragraph{时间序列包：xts}
\begin{verbatim}
library(xts)
data("sample_matrix")
sample.xts <- as.xts(sample_matrix)

sample.xts['2007']  # all of 2007
sample.xts['2007-03/']  # March 2007 to the end of the data set
sample.xts['2007-03/2007']  # March 2007 to the end of 2007
sample.xts['/'] # the whole data set
sample.xts['/2007'] # the beginning of the data through 2007
sample.xts['2007-01-03'] # just the 3rd of January 2007
\end{verbatim}


\subsection{缺失值的表格与图形探索}
\subsubsection{表格}
\begin{verbatim}
library(mice)
data(sleep,package='VIM')
md.pattern(sleep)
\end{verbatim}

\begin{figure}[H]
\includegraphics[scale=.8]{MissData0.png}
\end{figure}

缺失值有$ 42\times 0+2\times 1+3\times 1\cdots+1\times 3=38 $。


\subsubsection{图形}
\begin{lstlisting}
library(VIM)
aggr(sleep,prop=T,numbers=T)
\end{lstlisting}

\lstinline|prop,numbers|分别表示是否按比例显示，是否显示标签。

\includegraphics[scale=0.5]{MissData.png}
\subsubsection{检查面板数据是否平衡}
\begin{verbatim}
library(plm)
table(Grunfeld[,1:2]) # id 和 year 列
table(Grunfeld[-c(1,8),1:2])
\end{verbatim}

\section{画图：\texttt{ggplot2}}\label{chpDraw}
\texttt{ggplot2}其思路在于图层添加。（《R数据可视化手册》）
\subsection{例子1：各种标记的添加}
\begin{figure}[H]
\includegraphics[scale=0.6]{DrawDID.png}
\caption{共同趋势假设}\label{BR_1}
\end{figure}
相关代码见\verb|DrawDID.R|

\subsection{例子2：映射}
\begin{figure}[H]
	\includegraphics[scale=0.6]{scatter.png}
\end{figure}

\subsection{例子3：地图}
REmap包通过调用百度echarts的功能进行地图的绘制，绘制出的地图为html格式并具有交互功能。REmap包用来绘制地图的函数有四个：remap()，remapB()，remapC()和remapH()。四个函数绘制四种不同类型的地图，对导入的数据框的格式要求也不同。这里我们重点讲述使用remapC()函数，它主要是根据某一个指标，在各省份中用颜色标识填充，这在科研中应该是使用的频率较高的。

该包不是官网的标准包。因此安装稍费周折：
\begin{verbatim}
library(devtools)  
install_github('lchiffon/REmap',force = T) 
\end{verbatim}

看一个例子\verb|DrawMap.R|：
\begin{verbatim}
library(REmap)
province<-mapNames("china")    ##获取各省份名称
value<-runif(max=100,min=10,length(province))
popdata<-data.frame(province,value)

pop<-remapC(data = popdata, #分配各省的数值
theme = get_theme("Sky"), # 天蓝色主题
maptype = "china", # 可以是china, world, 或者中国省份的名字（参见popdata[,1]）
title = "中国各省市人口-随机")

plot(pop)
\end{verbatim}

在打开的网页中，右边有按钮可以输出图形。

\section{函数的编写}
\subsection{控制流}
\paragraph{循环结构}：{\color{RoyalBlue}\verb|for; while|}
\begin{verbatim}
for (i in 1:10){
   print('Hello')
}

i <- 10
while (i > 0){
  print('Hello')
  i <- i-1
}
\end{verbatim}
\paragraph{条件执行}
\begin{verbatim}
grade <- 'OK'
if (is.character(grade)) {
  grade <- as.factor(grade)
} else {
  print('grade is not a character')
}

if (is.character(grade)) print('OK') else print('Not') # 简化写法
\end{verbatim}

\paragraph{常用语句}
\begin{verbatim}
next()
break()
\end{verbatim}
\subsection{优化算法: \texttt{Rdonlp2}}
写一个MA(1)模型的条件极大似然估计，对于一个MA(1)过程，
\[ y_t=\mu+\varepsilon_t+\theta \varepsilon_{t-1} \]

其中$ \varepsilon_t\sim i.i.d.\; \mathcal{N}(0,\sigma^2) $。对于该模型，我们只有观测值$ y_t $，需要估计参数向量$\bm{\theta}= (\mu,\theta,\sigma^2) $。

要做极大似然估计，首先要写出$ y_t $的条件概率密度函数，这个比较简单，
\begin{align*}
y_t|\varepsilon_{t-1}\sim \mathcal{N}(\mu+\theta\varepsilon_{t-1},\sigma^2) \Longrightarrow f_{y_t|\varepsilon_{t-1}}(y_t|\varepsilon_{t-1};\bm{\theta}) & =\frac{1}{\sqrt{2\pi \sigma^2}}\exp\left[\frac{-(y_t-\mu-\theta\varepsilon_{t-1})^2}{2\sigma^2}\right]\\
& = \frac{1}{\sqrt{2\pi \sigma^2}}\exp\left[\frac{-\varepsilon_t^2}{2\sigma^2}\right]
\end{align*}

然后写出所有观测的联合似然函数，此处即为连乘，
\begin{align*}
\mathcal{L}&= f_{y_T,y_{T-1},\cdots,y_1|\varepsilon_0=0}(y_T,y_{T-1},\cdots,y_1|\varepsilon_0=0;\bm{\theta})\\
& = f_{y_1|\varepsilon_0=0;\bm{\theta}}(y_1|\varepsilon_0=0;\bm{\theta})\sum_{t=2}^{T}f_{y_T|y_{T-1},\cdots,y_1,\varepsilon_0=0}(y_T,y_{T-1},\cdots,y_1,\varepsilon_0=0;\bm{\theta})
\end{align*}

条件对数似然函数为，
\begin{align*}
\ln\mathcal{L}&=\ln f_{y_T,y_{T-1},\cdots,y_1|\varepsilon_0=0}(y_T,y_{T-1},\cdots,y_1|\varepsilon_0=0;\bm{\theta})\\
& = -\frac{T}{2}\ln 2\pi-\frac{T}{2}\ln \sigma^2-\sum_{t=1}^T\frac{\varepsilon_t^2}{2\sigma^2}\\
& = -\frac{T}{2}\ln 2\pi-\frac{T}{2}\ln \sigma^2-\sum_{t=1}^T\frac{(y_t-\mu-\theta \varepsilon_{t-1})^2}{2\sigma^2}
\end{align*}

我们的代码就是要基于观测值$ \{y_1,y_2,\cdots,y_T\} $选择合适的$ \bm{\theta}= (\mu,\theta,\sigma^2) $最大化上式。写这个代码要注意的就是全局变量的使用。

代码见\verb|MLE.R|
